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Abstract 

This paper deals with the question of how to most effectively con- 
duct experiments in Partially Observed Markov Decision Processes so 
as to provide data that is most informative about a parameter of in- 
terest. Methods from Markov decision processes, especially dynamic 
programming, are introduced and then used in an algorithm to maxi- 
mize a relevant Fisher Information. The algorithm is then applied to 
two POMDP examples. The methods developed can also be applied 
to stochastic dynamical systems, by suitable discretization, and we 
consequently show what control policies look like in the Morris-Lecar 
Neuron model, and simulation results are presented. We discuss how 
parameter dependence within these methods can be dealt with by the 
use of priors, and develop tools to update control policies online. This 
is demonstrated in another stochastic dynamical system describing 
growth dynamics of DNA template in a PCR model. 

1 Introduction 

Hidden Markov Models have proven their usefulness across a wide variety of 
applications. In many of these applications, the user or the experimenter will 
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have some way of influencing the transitions of the underlying Markov Chain, 
as in Markov Decision Processes, and such a process is called a Partially 
Observed Markov Decision Process, see Monahan [3]. If we assume that the 
transition probability matrix is governed by some unknown parameters, an 
important problem is to understand how the process can be influenced to get 
data that is most informative about the parameters. We can think of this 
as experimental design for Partially Observed Markov Decision Processes 
(POMDP). 

We consider a POMDP (x t ,yt,u t )t=i,...,T- In this setting x t is an unob- 
served Markov Chain, where the transition probabilities depend in a para- 
metric way on what control u t is chosen at time t and an unknown parameter 
9. The process y t is observed and depends on which state x t is in. 

Our goal is to find ways to use the controls Ut to improve parameter es- 
timates of 9. Since maximum likelihood estimates for 9 are asymptotically 
efficient, our general strategy is to use the controls to try to minimize the 
sample variance of maximum likelihood estimates of 9. This will be achieved 
by maximizing a Fisher Information for 9. The controls are calculated us- 
ing dynamic programming, a popular maximization algorithm from Markov 
Decision Processes which outputs an adaptive control policy, i.e. the control 
chosen at time t is based on observations up to time t. 

Lin et al. [2] proposed maximizing the Fisher Information that corre- 
sponds to direct observations of the entire process x t , labeled the Full Infor- 
mation Fisher Information (FOFI), and using a filter to compute Xt if it is 
not observed directly. This paper makes use of the POMPD structure and 
proposes maximizing a Fisher Information that is based on the observations 
y t , labeled the Partial Observation Fisher Information (POFI). The control 
policies based on these two Fisher Informations, FOFI and POFI, are com- 
pared, and we illustrate a setting in which the two control policies are quite 
different. 

The methods we use to calculate controls for maximizing Fisher Infor- 
mation will depend on the unknown parameter 9. We illustrate how this 
problem can partially be overcome by assuming a prior for 9 to calculate a 
control policy before running the experiment. Additionally we describe how, 
using data acquired as the experiment progresses, a posterior for 9 can be 
used to calculate a more precise control policy. That is, parameter informa- 
tion from observations acquired at some time t can be used to improve the 
policy used in what is left of the experiment. These methods will be based 
on the Value Iteration Algorithm (VIA), which is closely related to dynamic 
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programming. 

The methods developed have application value beyond Partially Observed 
Markov Decision Processes. Lin et al. [2] considered stochastic systems of 
the form 

dx = f(x, 9, u{t))dt + S 1/2 dW 

where 9 is the parameter of interest, to be estimated, u(t) is a control that 
can be chosen by the user, x is the vector of state variables, f is a vector 
valued function, W a Wiener process, and additionally x(t) is only observed 
partially or noisily. By discretizing time, state and observation spaces the 
process can be approximated by a POMDP, allowing us to use the methods 
developed in this paper to devise a control policy that maximizes information 
about the parameter 9. 

In order to illustrate these methods we present four examples, two par- 
tially observed Markov decision processes, and two continuous stochastic 
systems, that need to be discretized before our methods are applied. In the 
first three examples we use the unknown 9 to calculate our controls, but in 
the last example we compare using a prior for 9 to calculate a control policy, 
and using VIA to calculate a control policy that is updated online. 

In the first example we hypothesize about the kind of systems in which we 
will observe large improvement in parameter estimation by using the POFI 
control policy over the FOFI policy. Following a discussion we construct a 
mock Partially Observed Markov Decision Process, in which this improve- 
ment is shown using a simulation study. 

To illustrate the real- world applicability of design in discrete POMDP 's 
we consider a realistic POMDP from experimental economics. The model will 
consist of a simple adversarial game similar to the "rock - paper - scissor" 
game where one player tries to play in such a way that maximizes information 
about the other players' strategy. 

The third example is a stochastic version of the Morris-Lecar Neuron 
model, a dynamical system which models voltage in a single neural cell. 
This model is two dimensional, but only one dimension is observed. The 
model has multiple parameters and we investigate how the POFI and FOFI 
control policies perform in estimating them. 

Fourthly we consider an example from biology, a Polymerase chain reac- 
tion (PCR) experiment where DNA template is grown in liquid substrate. 
The population dynamics are modeled in a dynamical system with stochastic 
errors, and the aim is to estimate the half-saturation constant, a parameter 
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which controls the saturation of the template. Here we compare using a prior 
for 9 and using VIA to calculate a control policy. 

In Section 2 we describe the POMDP setup, and then introduce dynamic 
programming, an algorithm from Markov Decision Processes. This algorithm 
is then applied to evaluate control policies that maximize both the Partial 
Observation Fisher Information and the Full Observation Fisher Information. 

In Section 3 we show how the control policies can be applied to POMDP 
examples and then, in Section 4, describe how stochastic dynamical systems 
can be transformed to POMDPs by using discretization methods, and give an 
example of such a system. In Section 5 we discuss the problem of parameter 
dependence when running the dynamic programs, and how these issues can 
be dealt with. This is illustrated in the last example, and a discussion follows. 

2 Basic Setup 

We consider a Markov decision process (x t ,u t )t=i,...,T- In this setting x t is 
a Markov chain, but the transition probabilities at time t depend on the 
control u t chosen at that time. We assume a finite state space S for the state 
process x t and that the controls available belong to some finite set A. The 
transition probability matrix P which depends on a parameter 9, and the 
control u r chosen is given as 

P[j(0) = p(x t+1 = x l \x t = x 3 ,u t = u r ,9) 

where x\ x j G S and u r G A. In addition to this we assume that the process 
Xt is latent and we only observe the related observations yt whose relation to 
the x t process is described through the matrix E with 

Eij =p(y t = y l \x t = x 3 ) 

This makes the system a Partially Observed Markov Decision Process. It has 
a finite horizon T in which we observe Yt = Hi ■ ■ ■ Ut- The objective is to 
use the controls u t to maximize the information we get about the parameter 
9 through the observed process yx- The parameter estimation is done using 
maximum likelihood and it is therefore natural to try to maximize the Fisher 
Information of our observed process 

FI{9) = E[-l"{9)\ = E 
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We label this as the Partial Observation Fisher Information (POFI). When 
we consider continuous time dynamical systems the observation spaces will 
be continuous, but we will use this discretized Fisher Information as an ap- 
proximation to the actual Fisher Information of the observations. 

2.1 A Dynamic Program 

In order to maximize the Fisher Information we employ the techniques of 
stochastic dynamic programming. We assume a Markov Chain (xt)t=i,...,T, 
whose transitions are influenced by controls {ut)t=i,...,T-i, an associated re- 
ward function C t (x t ,u t ) and our objective is to maximize the total expected 
reward E[J2j =1 C t (x t ,u t )} by use of the controls. The essence of dynamic 
programming is that by starting at time T — 1 and working backwards, we 
can compute an optimal policy to map Xt to Ut that accounts for the choices 
of u t that we will make in the future. Starting at time T — 1 we calculate 
the optimal reward to go: 

V T -i(x T -i) = m.ax{E XT [C T -i(xT-i,UT-i)\x T -i,UT-i\} 

UT-l 

and the associated control 

«t-i(«t-i) = argmax{S a!r [C7 r _i(arr-i,ttT-i)|a;T-i,«T-i]} 

UT-l 

Both are calculated for every state xt-i- Then, working backwards through 
t = T — 2, . . . , 1 we solve the optimality equations 

V t (x t ) = max{E x [C t (x t ,u t ) + V t+1 (x t+1 )\x t ,u t }} 

Ut 

and get the associated control 

u* t (x t ) = axgmax{E Xt+1 [C t (xt,u t ) + V t+1 (x t+1 )\x t ,u t ]} 

u t 

for every state Xt- This will give us a policy of what control to use at a 
certain state xt at a certain time t. The use of these controls will maximize 
the total expected reward E[J2 t Ct(x t ,Ut)]. We refer to Puterman [5] for a 
detailed description of dynamic programming. 

With this in mind, we write y t = (yt, ■ ■ ■ , Vi) and u t = (u t , ■ ■ ■ , U\) and set 
C^y^Ut,^) = -J^j logp(?/ m |yt, u t , 6) and we try to maximize the Fisher 
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Information FI{9) = J2t E[C t {y t , u t , 9)\y t , u t , 9}. Note that in this instance 
the cost function depends on the entire history of observations y t and controls 
u t . 

Define the Fisher Information To Go (FITG) to be 



V5(y t ,ut_i,0) = max{^ w+1 [C t (y t ,Ut,0) + V t+1 (y t+1 ,u t ,9))\y t ,u t ,9}} 
logp(y t +i|y t ,u t ,0)) + V t+1 (y t+1 ,u t ,9)\y t ,u t ,9 



max <J E yt+1 



u t 

d 2 



89 
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Pseudocode for this Dynamic program can now be given as follows 
V T = 

for t = (T—l)—>l do 

V y t and u t _i calculate and store 

V t (yt,ut-i,6) = max ut {E yt+1 [C t (y t , u t , 9) +V t+1 \y t ,u t ,9]} 
u*(y t ,u t -i,9) = argmax Ut {E yt+1 [C t (y t , u t , 9) + VJ+i|y t , u t , 9] } 
end for 

This process will provide us with a best control, a best action at some 
time t, for every past history of states and controls and represents, if it can 
be calculated, a gold standard for maximizing the Partial Observation Fisher 
Information. 

A problem here is that in the first step of the dynamic program (t — T— 1) 
we would have to calculate the FITG for L T ~ l l T ~ 2 many combinations of y t 
and u t _i, where L is the dimension of the y t state space and / the number of 
controls we can choose from at each time. This will be formidable for even 
modest dimensions. We therefore approximate the process by using only the 
values for the past two time points within the Fisher Information. 

P{yt+i |yt, u t , 9) w p(y t +i\y t , u t , y t -i, iit-i, 9) 

and concentrate on maximizing an approximation of the Partial Observation 
Fisher Information 



FI = E 

± ± approx / J- 1 ] 



Vt+l 

t=0 



& 

— logp(y t+1 \y t , u t , y t -i,ut-i, 9)\y t , u t , y t -i,u t -i 



It is possible to consider longer sequences of past observations if needed, but 
at the expense of computational time and storage. 
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With this approximation, we set the Fisher reward to be 

d 2 

C(y t ,ut,yt~i,u t -i,9) = -— logp(y t+l \y t , u t , y t _i, u t -i, 9) 
and the Fisher Information To Go to be 

V t {y u y t -uu t -x,0) = max [E y [C + V t+1 \y u u u y t ^u t ^9)\] 

As obtained above the pseudocode for this dynamic program is: 
V T = 

for £ = (T — 1) — >• 1 do 

V y t ,yt-i,u t -i and calculate and store 

Vt(yt,yt~i,u t -i,9) =max Uf {E yt+1 [C + V t+1 \y t , u t , y t -i, u t -i, 9] } 
u$(yt,yt-uu t -i,9) = argmax Ut {E yt+1 [C + V t+ i\y t , u t , y t -i, «t_i, 9)} 



Explicit expressions for C(y t , u t , y t ^ x ,u t -i, 9) and p{y t +i\y t , u t , y t -i,«t-i, &) 
are given in the Appendix. 

2.2 FOFI Dynamic Program 

An alternative method to choose controls was proposed by Lin et al [2]. Here 
we construct an optimal control policy for the Fisher Information that would 
apply if (x t ) were observed directly. 



This is labeled the Full Observation Fisher Information (FOFI). As noted 
before, when considering continuous time stochastic systems, the state space 
will be continuous, but we will use this Fisher Information as an approxi- 
mation to the continuous state Fisher Information. The advantage of using 
FOFI over POFI is that when running the dynamic program the Markov 
property of Xt allows us to only consider a maximization over the state space 
of x t but not past values x t -i,x t -2, 

Here the reward function is C(x t ,u t ,9) = —-^2logp(x t+ i\x t ,u t ,9). The 
pseudocode for this dynamic program is: 



end for 




t=i 
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V T = 

for t = (T - 1) 1 do 

V x t and calculate and store 

V t (x t ,9) = max Ut {E Xt+1 [C(x t ,u t ,9) + V t+1 (x t+1 ,9))\x t ,u t ,9]} 
u*(x t ,9) = argmax Ut {E Xt+1 [C(x t , u t , 9) + V t+1 (x t+1 ,9)\x t ,u t ,9]} 
end for 

However, the problem with FOFI is that when the actual experiment is 
run we don't observe x t . Instead we have to use the observed values to get a 
probability distribution (a filter) on x t , p(xt\yt, Ut-i, 9) and use the control 
associated with the state that has the highest probability. 



3 Discrete Examples 
3.1 6 state example 

While the FOFI strategy has been shown to be effective in Lin et al. [2] it is 
possible to define systems in which it is not optimal and may in fact be worse 
then just using fixed or random controls. In some systems, certain parts of 
the state space give more information about a parameter than others, and it 
would make sense for the controls to push the process towards these areas. It 
does matter, however, how well these parts of the space are observed; if they 
are observed with a great deal of noise, or not at all, it might be suboptimal 
to push the process there. It is in cases like this that POFI often does better 
than FOFI. In this example, we demonstrate a system where FOFI and POFI 
choose very different controls, and using a simulation study, we show that 
the controls chosen by POFI produce less variable parameter estimates. 

Consider a discreet time Markov chain x t with state space S x = {1, 2, 3} 
and a transition probability matrix 



P 



!' 



+ ~ 



1 
2 

E 
2 

1 p u 

2 4 4 



4 _ « 
4 

.15 

.45 + f 



where the parameter of interest is p G [0, .5] and the control is u G {—1, 1}. 
For Xt = 1 or x t = 3, choosing the control u = 1 will increase the probability 
of the Markov chain staying in its current state while choosing u = — 1 will 
increase the probability of it leaving its state. 
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Ut 


Vt 


Vt-i 


«t-i 


1 


1 


1 


1 


-1 


2 


1 


1 


-1 


1 


2 


1 


1 


2 


2 


1 


-1 


1 


1 


-1 


1 


2 


1 


-1 


1 


1 


2 


-1 


-1 


2 


2 


-1 



Table 1: Long run control policy that results from using POFI in the 6 state 
example. The first column describes which control to use for a given history 
(y t ,yt-x,Ut-i) of observations and control. 

Now assume this process isn't observed directly but through a related 
process y t with state space S y = {1, 2} whose transition probabilities depend 
on which state Xt is in. We denote the transition probability matrices with 
(Pk){i,j} = p(yt+i = j\vt =i,x t = k) given by 



If xt were observed we would get information about the parameter p when xt 
leaves state 1 and from y t when xt = 3. The idea here is that since the FOFI 
controls assume the whole state space is observed they might encourage x t 
to be in state 1, while the POFI controls that take into account what is 
actually observed might choose the controls more intelligently. Indeed when 
calculating the controls according to FOFI the long run control is to "leave 
one's state" if xt — 3 and "stay in one's state" if Xt — 1. It is harder to 
predict what kind of controls result from using POFI, but the control policy 
is given in Table [T] 

To illustrate this difference, a simulation study was carried out to test 
which method performed better: The process Xt was run for 1000 steps with 
p = .37, using controls chosen by POFI and again using those chosen by 
FOFI. Additionally we ran a simulation of the same length, but where the 
control was chosen randomly, with u = —1 and u = 1 having equal proba- 
bility. Then the parameter p was estimated using an EM algorithm. This 



.5 .5 
.5 .5 



I- 2 

2 

2 



P 
2 

1 - 2 
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bias st. dev. MSE 


FOFI 
POFI 

Random 


.0109 .1098 .0121 
.0094 .0867 .0076 
.0059 .0941 .0089 



Table 2: Simulation results for the 6 state example. We see that the controls 
chosen by POFI make for more accurate estimates of p. The FOFI policy 
does worse than a random policy. 





Left 


Right 


Left 


2,0 


0,1 


Right 


1,2 


1,1 



Table 3: Rewards in the Gamble Safe game. The first number is the reward 
for the Row player and the second number the reward for the Column player, 
given a certain outcome. 



was done 500 times to get an empirical distribution for the estimates of p. 
The results are given in Table [2] Estimates of p using a POFI policy had the 
lowest MSE and variability, but estimates based on the random policy had 
lower bias. Estimates using a FOFI policy were comparatively worse. 

3.2 Adversarial game - A POMDP example 

The following example describes how the methods above can be used in the 
context of experimental economics. In the problem described below we are 
interested in eliciting the Gamble-Safe game, described in Sachat et al. [6] 

This system is a game with two players: a Row player and a Column 
player. They repeatedly play a game where both simultaneously choose either 
left or right, and they get rewards depending on the outcome, according to 
Table [3] the Row player would for example get 2 and the Column player 
if both chose left. We follow [6] and assume that at any given play the 
Column player follows one of two strategies: the Nash-equilibrium strategy of 
choosing either left or right with 50% probability or the Gamble-safe strategy, 
where they only choose right. The player will pick either strategy based on 
a multinomial logistic model, where the probabilities depend on the last 
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two plays of the Row player, and the last strategy chosen by the Column 
player. This results in a Partially Observed Markov Decision Process with 
the strategy employed being a hidden state giving rise to observed plays. 

Let S t denote the strategy chosen by the Column player at time t, U t 
denote the action played by the Row player at time t. Let S t — — 1 if the 
Nash-equilibrium is chosen, S t — 1 if the Gamble-safe strategy is chosen. 
Also let Ut — 1 if the Row player plays right, U = — 1 if he plays left. 
Similarly Y t will denote the plays of the Column player. The strategy St+i 
chosen at time t + 1 will then be chosen according to 

= -1) = ^ and P(S t+1 = 1) = ^ 

where we let 

x = l.2U t + U t -! + 6 S t 

Our objective is now as follows: assume the Row player knows what kind 
of model the Column player uses to pick strategies, is interested in estimating 
9, and doesn't really care as much about winning. How should they play in 
order to estimate 9 with maximal precision? 

To cast this into our usual setting we think of S t being the unobserved 
underlying Markov Chain, U t as the control and Y t as the observed process. 
Since the transition probabilities from St depend on Ut-i (a part of the history 
at time t — 1) we augment the state space to include U t -\, i.e. R t = (S t , U t -\) 
will be our underlying Markov Chain. At this point we could run the dynamic 
programs for both FOFI and POFI, but controls calculated that way will 
depend deterministically on the plays of the Column player. Seeing that 
realistically deterministic plays can often easily be countered in adversarial 
games, it is better to follow a strategy that includes some randomness in the 
plays. So we let W t G {—1, 1} be the strategy of the Row player in such a 
way that 



% = \ W ' P ' o ) *Wi = l> and % = \ W ' P ' '1 ) ifWi = -l 
U t = -1 w.p. .2 J ' U t = -1 w.p. .8 J 

These kind of changes are easily incorporated in the dynamic program for 
both POFI and FOFI, by adding an expectation over W t at every step t. 
We give the pseudocode for the FOFI dynamic program as an example. Set 

C(r t ,u t ,9) = -J^logp(r t+ i|r t ,ut,0). 



11 





Ut-1 = 1 Ut-2 = -1 

Y t .! = i Ft_i = -i y t _i = i y t _i = -i 


u - 1 Yt = 1 
^ - 1 y t = -l 

u - a Yt = l 

U ^-~ l Yt = -1 


111-1 
-1 -1 -1 -1 
-1 -1 -1 -1 
1111 



Table 4: Long term POFI control policy: What U t to choose given different 
combinations of Y t , Y t -i, U t -i, Ut-2- 



V T = 

for t = (T - 1) 1 do 

V r t and calculate and store 

V t (r t , 9) = max wt {E Ut [E n+1 [C + V t+1 (r t+u 6))\r t , 6} \w t ] } 
w* t {r t , 9) = argmax^ {E Ut [E n+1 [C + V t+1 {r t+1 , 6))\r t ,6] \w t ] } 
end for 

We set 9 = .7 and calculated the FOFI and POFI policies. The long run 
FOFI policy was to alternate between W t = 1 and W t = —1 at every step 
(no matter the state of Rt). Thus if Wt — 1 at time t, we choose Ut = 1 with 
probability .8 and then at time t + 1 we will set Wt+i — — 1. The long run 
POFI policy depended on F t , Yj_i, Z7t_i, C/t_2 in a rather complicated way, see 
Table SI 

To compare the two policies we ran a simulation study with T = 500, 
and 500 simulations for both the POFI and the FOFI policy. Another 500 
simulations were run where the plays (control) where chosen randomly. The 
parameter 9 was estimated using an EM algorithm. The results of this es- 
timation under each policy are given in Table [5] where the POFI controls 
produce least variance and the most accurate estimates. 

4 Discretization methods 

In order to apply the methods described above to dynamical systems, we 
need to approximate them by a suitable Partially Observed Markov Decision 
Process. We achieve this by discretizing time, state and observation spaces. 
In this paper, the continuous stochastic dynamical systems considered are of 
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bias st. dev. MSE 


FOFI 
POFI 

Random 


0.004 0.261 0.068 
0.016 0.231 0.054 
0.013 0.252 0.064 



Table 5: Simulation results for the adversarial game. The POFI control 
policy seems to offer the best controls for optimizing the estimation of 9. 
FOFI performs similarly to a random policy. 



the form 

d* = t(-K,6,u(t))dt + Y^ 2 (W 

where 9 is the parameter of interest, to be estimated, u(t) is a control that 
can be chosen by user, x is the vector of state variables, f is a vector valued 
function and W a Wiener process. The dynamical system is approximated 
on a fine grid of times i t — t5, t — 1, . . . , T and we obtain a discrete-time 
model 

x t+1 = x t + <Jf(xt, 9, u t ) + VSe lt 

where e u ~ iV(0, Si) are independent normal random variables. We assume 
the underlying state variables x t are only observed partially or noisily. 



Yt = g( x *) + e 2< 



where e 2t ~ iV(0, £|). 

In order to approximate this as a Markov Chain, the state space is dis- 
cretized in each dimension and the model is then thought of as moving be- 
tween the different boxes. The probability of moving from box to box is 
approximated using the normal pdf at the midpoints of the boxes. In the 
examples in this paper, only equidistant discretization is considered, but this 
restriction can be readily removed. If we label the two midpoints as %\ and 
ii and the area of the second box as area x this probability is given as 

p(x t+ i = i 2 \x t = ii,ut,9) 
1 

(27r) fc / 2 det(S 1 ) 1 /2^^ 2 ( 

where k is the dimension of x. The probabilities are then normalized to 
make sure they sum to 1. If the controls u t can be chosen on a continuous 



exp --(ii + 5f(ii,0,u t )) (h + Sf(i 1: 9 : u t )) )-area x 
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scale then this scale has to be discretized as well. (x t ,u t ) is then a Markov 
Decision Process, and one can run the FOFI dynamic program. 

For the POFI dynamic program the observation space needs to be dis- 
cretized as well. The probability of what observation box is observed depends 
on in which box the underlying Markov Chain is in. If we label the midpoint 
of the underlying Markov chain midpoint as % and the midpoint of the ob- 
served process box midpoint as j, and the area of the latter box as area y this 
probability is given as 



These probabilities are also normalized to sum to 1. The process (x t ,yt,u t ) 
is now a Partially Observed Markov Decision Process and one can run the 
POFI dynamic program. 

4.1 Morris Lecar Model 

The Morris Lecar Model [7j describes oscillatory electric behavior in a single 
neural cell, as regulated by flow of Potassium and Calcium ions across the 
cell membrane. These models are defined in terms of state variables vt and nt 
representing the voltage across the membrane and the flux of the Potassium 
channel respectively. 



C m v t = I t -gr (v t - Ei) - g K ■ n t ■ (v t - E K ) - g Ca ■ m 00 (^) • (v t - E Ca ) 



p(yt = j\x t = i) 



1 



-gO' - 9(i)) s 2 U - 9{i)) j -area. 



(2vr) fc / 2 det(S 2 ) 1 /2 



cxp 



y 



(1) 



n t = -(f) ■ (n t - n^ivt)) / r n {v t ) 



(2) 



where 




T n{ v ) — sech((i; — t>3)/(2t> 4 )) 
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The voltage between cells depends on Potassium and Calcium concentrations, 
and on the amount of leakage. The further these factors are away from 
their equilibriums Ei, Ek, Ec a the greater the rate of change in voltage. The 
multiplicative value n t which has to do with Potassium levels, is modeled 
through the second differential equation in which n t is driven towards a 
voltage-dependent equilibrium level defined by n^iyt) but converges to this 
at a much slower rate then the dynamics of Vt- The neuron is stimulated by 
an external applied current, It (our control), and Vt is measured. Our goal 
is to maximize information about the parameters C m ,gc a and 0, considered 
separately. 

We consider a stochastic version of this model 

C m dvt = 

It- 9r (v t - Ei) - g K -n t - (v t - E K ) - g Ca ■ m oc (^) • (v t - E Ca ) + a ■ dW x 

dn t = • (raoo(i>i) - n t )/r n (v t ) + a ■ dW 2 
and the first step is to discretize these equations with respect to time: 

v t (t + dt) = 

v(t)+dt(I t -g l iv t -E l )-g K m t iv t -E K )-gca-rn^(v t )iv t -E C a))/C m +aVdi-e 1 

n t (t + dt) = n(t) + dt((f) ■ (n^Vt) - n t )(t)/T n (v t )) + a\fdt ■ e 2 

where e±, e 2 ~ N(0, 1). 

We discretized v t onto the range [—75, 45] and n t onto [0, 1], after running 
a few trial versions of the model. Both ranges where discretized into 25 
intervals. Only v t is measured and it is measured noisily, 

y t = v t + e t 

where e t ~ iV(0, 1). The observation space was discretized to the same range 
as Vt but into 20 intervals. These approximations give rise to a Partially 
Observed Markov Decision Process to which our methods can be applied. 
Both FOFI and POFI controls were calculated for this experiment, for the 
parameters C m , gc a and 0, considered separately. The values for the param- 
eters were set to be C m = 20, g Ca = 4.4, g t = 2.0, E k = -84.0, E t = -60, 
E Ca = 120.0, = .04, Vl = -1.2, v 2 = 18.0, v 3 = 2.0, v 4 = 30.0, a = a = 1 
and dt = 1. The controls range was set to be [—1.5, 6.0] and discretized to the 
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parameter 




bias 


j_ i 

st. dev. 


iv inn 

MSE 




r Ur 1 


A OQ A 

AZ64: 


O A 700 

2 A ( 22 


p. on 1 o 


r 


POFI 


.4129 


2.4068 


5.9632 




Fixed 


.9098 


3.4240 


12.551 




FOFI 


.0613 


.3671 


.1385 


9Ca 


POFI 


.0158 


.3706 


.1376 




Fixed 


.0249 


.6193 


.3841 




FOFI 


.00485 


.01085 


.00014 




POFI 


.00257 


.01037 


.00011 




Fixed 


.01357 


.02643 


.00088 



Table 6: Simulation results for the Morris-Lecar model, consider the pa- 
rameters C m ,gca,4> separately. We see that the POFI and FOFI policies 
outperform the fixed policy It = 1.5 in all cases, and the POFI policy seems 
to perform slightly better than the FOFI policy for the three parameters 
considered. 



set I t G { — 1.5,0.0,1.5,3.0,4.5,6.0}. An example of what the control policy 
looks like is given in Figure [T] A simulation study was run for each of the 
three parameters; the system was simulated within the discretized Markov 
Chain framework with 100 time steps and all schemes had 100 simulations. 
The parameter in question was estimated for each simulation using an EM 
algorithm. As a baseline comparison we also ran a simulation study using 
a fixed control (I t = 1.5). The results are given in Table [6] The difference 
between POFI and FOFI turns out to be not very dramatic, likely due to 
the observations providing a great deal of information about the underlying 
state variables, which is when FOFI performs well. 

5 Parameter dependence of dynamic program 

Until now we have calculated the dynamic program assuming knowledge of 
the parameter 9, the very thing we wish to estimate with maximal precision. 
Since the dynamic programs we have considered are run before the exper- 
iment is started we generally won't have data to estimate 9. Additionally, 
for the FOFI simulations we have used 9 directly to estimate Xt within the 
filter to get the appropriate control, but this will not be possible in practice. 
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next control control next control control 




voltage v y_t 



(a) FOFI control (b) POFI control 

Figure 1: Long term controls of FOFI and POFI for the parameter gc a - The 
FOFI plot gives the control to use, given a certain position in state space 
(v t ,n t ). The POFI control to use will depend on the last two observations 
and the last control, but fixing the last control as, for example, J t _i = 6 one 
can plot which control to use given combinations of last two observations. 
In both the FOFI and POFI case, the control chosen is mostly either at the 
extreme high or extreme low end of the range of controls 
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There few ways of dealing with this. 

Assuming some prior information one can use a prior for 9 to run the dy- 
namic program. To do this, we add one more expectation for 9 at every time 
step t, and then maximize the expected Fisher Information to get the best 
control. In the POFI case this maximizes the expected Partial Observation 
Fisher Information: 



T-l 



E e [FI{9)] = Y J Ee 



t=i 



E 



yt+i 



(J 2 



89' 



logp(2/i + i|yt,ut,0)|yt,u t ,0 



or more practically, the expected approximated Partial Observation Fisher 
Information: 



Eg[FI, 



approx 



T-l 



t=l 



E 



yt+i 



2 



d9 2 



logp(y t+ i \y t , u t , y t _i, u t -i,0) \y t , u t , y t -i,u t ^ 9 



In the FOFI case this would maximize the expected Full Observation Fisher 
Information: 



T-l 



E e [FI] = Y J Ee 



t=i 



Ext+i 



d 2 



:p(x t+1 \x t ,Ut,9)\x t ,u t ,9 



This strategy was employed in Lin et al. [2]. 

The rather obvious deficiency here, for both POFI and FOFI, is that as 
the experiment runs, we get observations that can be used to improve our 
prior for 9, and could be used to get better controls, if we could brake the 
experiment and rerun the dynamic program. 



5.1 Online updating 

In some systems the time spent in each state is very short, too short to 
perform many calculations, making it valuable to have a "look-up table" of 
controls. Here the POFI controls have an advantage over the FOFI controls, 
in the sense that they are of the "look-up" kind, while FOFI requires esti- 
mation of the underlying x t process, before the control can be looked up, 
usually done using a filter. 

In other systems, there is time to do some calculations between transi- 
tions. Note, for example, that at time t we have observed yi,...,y t and 
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this will allow us to calculate a posterior distribution 7r(0|y t , u t _i) for our 
parameter of interest. This posterior could then be used to run the dynamic 
program again, as described above, from time T — 1 to time t. This can be 
quite time consuming if done at each time step t, so we propose a method 
that relies on the Value Iteration Algorithm. 



5.1.1 Value Iteration Algorithm 

A popular algorithm from the theory of Markov Decision processes is the 
Value Iteration Algorithm (VIA), see Puterman [5]. The theoretical moti- 
vation of the Value Iteration Algorithm is similar to the dynamic program 
described above, but here the objective is to maximize an expected total 
reward W\ that has a discounting factor A, where < A < 1, and the time 
horizon is assumed to be infinite; 



W 1 = E 



t=i 



and Wi is labeled as the expected total discounted reward. In Puterman [3] it 
is shown that an optimal control policy exists and it can be chosen in such a 
way that the control is time independent, i.e. depends only on the state of Xt, 
and not the time t. Moreover this optimal control can be approximated using 
the Value Iteration Algorithm described below. Our experimental setting 
will neither be discounted nor have an infinite time horizon, but another 
result from Markov Decision Processes is that the control that maximizes 
the expected total discounted reward also maximizes the expected average 
reward W2 where: 



Wo 



lim — E 

n— »oo fl 



'N+n 
t=N 



Xt, Ut) 



given that A is chosen above a certain (unknown) threshold, see [S], although 
choosing A too high will cause the algorithm to converge slowly. The standard 
pseudocode for VIA is 



Set v 1 = and n = 1 
while \ \v n — f n_1 || > £ do 
V x t and calculate and store 

v n+1 (x t , 6) = max u {E Xt+1 [C(x t ,u t , 6)) + A • v n (x t+ll 9))\x t , u t , 6} } 
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n=n+l 
end while 

u*(x t ,9) = argmax Ut {E Xt+1 [C t (x t ,u t ,9) + A • v n (x t+ i, 9)\x u u t , 9] } 

where we label v n as the value vector. Since the main operation of VIA is 
a contraction mapping, the algorithm runs until it converges to some fixed 
point, within some error. 

Our aim with VIA is to maximize what we in the POFI case label, the 
average Partial Observation Fisher Information 



or in the FOFI case, the average Full Observation Fisher Information 



which is a reasonable quantity to maximize in order to obtain a time-invariant 
policy. 

We propose running VIA at every time step t, but to use the posterior for 
9, 7r(0|y t , u t _i), which is conditioned on all the data observed so far, instead 
of using the prior for 9. This will give a control that maximizes the average 
Fisher Information, using all the parameter information that is available at 
time t. Instead of starting VIA at each time t with v 1 = 0, considerable time 
can be saved by using the last value vector v n from the previous run of VIA 
at time t — 1. This is because the posterior for 9 often doesn't change much 
between time steps, and the last v n from time t — 1 thus being relatively 
close to the fixed point at time t. We now describe the pseudocode for this 
experimental design, using POFI. 

Let v ™ denote the value vector at time t at the n'th iteration of the t'th 
VIA and let 7r(6>|y t , u t _i) denote the posterior for 9 given observations up till 
time t. Also, to ease notation, let z t = yt, yt-i, Ut-i- 

Set v Q x = and n = 
for t = 1 T do 

while | \v n — v n ~ l 1 1 > e do 
V z t and calculate and store 
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< +1 (zt) 



d 2 



max { V ^ logp(y t+ i |z t , «t, 0) + Aw t n (z t+ i))p(?/ t+ i |z t , w t , 0)7r(0|y t , u t _i) 



Hi 



8 yt+i 



n=n+l 
end while 

Set = v? 
Now let 
«t(zt) = 

{d 2 
zC Z-^ - ^ lo SP(^+i l z t, «t, 0) + Aw t n (z t+1 ))j9(y t+1 |z t , itt, 0)7r(0|y t , u t -i) 
6 yt+i 

Use control u t , and observe and then update the posterior for 0, 

(a \ \ P(y*+i|yt,ut,0)7r(0|yt,u t _i) 
7r(0y t+ i,u t J = ^ — -. ; ^-t^ r 

end for 

Updating FOFI policies online using VIA can be done in a similar way. 
In the next example we compare using updated policies. 



5.2 PCR model 

Polymerase chain reaction is a well established method to copy and multiply 
DNA. We are interested in modeling the growth dynamics of DNA template 
(x t ), for a fixed amount of substrate. The model we use is 

g(1 — Ut)xt r-r 

x t+ i = (1 - u t )x t + at - — — + Vat ■ e\ 

(b + (1 - u t )x t y 

where Ex ~ N(0,af). Here Xt is the amount of DNA template, a and b 
the parameters of the model and u t the control, the percentage of template 
removed at each time point. We are interested in estimating the parameter 
b, labeled the half-saturation constant. A good reference for PCR models is 
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We measure the amount of DNA template at each time point, but with 
an error. Our observations are 

y t = x t + e 2 where e 2 ~ N(0, of) 

and thus we have a dynamical system which when discretized becomes a 
Partially Observed Markov Decision Process. 

The range for xt was set to be [0, 15] and then discretized into 200 in- 
tervals, and y t was discretized to the same range, but only into 20 intervals. 
The parameter values were set to be a = 2.0, b = 4.2, a± — a% = 1, dt = 1 
and the possible values of the control u t G {0, .2, .4, .6, .8, 1}. 

The objective is to maximize both the Full Observation Fisher Infor- 
mation, and the Partial Observation Fisher Information, but here we more 
realistically replaced the correct parameters with priors. We conducted a 
simulation study using controls based on these priors for both POFI and 
FOFI, and then compared their performance to controls that are updated 
online using VIA, also both for POFI and FOFI. As a baseline comparison 
we also ran a simulation using a fixed control. 

The range for b was set to be b G [1.7,8.0] and then we discretized that 
interval into 10 points {1.7,2.4,3.1,3.8,4.5,5.2,5.9,6.6,7.3,8.0}. We then 
considered a uniform prior on these points and a prior that is somewhat 
inaccurate, and puts the weight .9 on the point 7.3 and gives the others 
equal weight. The discounting factor for VIA was set to be A = .9. 

Our simulation study had the time length T = 200 and there were 300 
simulations for each case. The parameter b was estimated using an EM 
algorithm. The simulation results are given in Table [7] 

We note that when we calculate the controls prior to the experiment (No 
online updating), both the POFI and FOFI controls are significantly better 
than using the fixed control, and POFI seems to do a bit better the FOFI. 
Interestingly, calculating the controls using the inaccurate prior does better 
then using the uniform prior, likely due to a reduction in variance. 

Accuracy increases in most cases when we allow for online updating using 
the VIA algorithm. Starting the VIA with an uniform prior does better than 
starting with the inaccurate one, which is probably due to the VIA having 
to spend more time fixing the prior. 

Additionally, in Figure [2j we see that using the previous final value vector 
as the starting value vector of VIA when going from time point t to t + 1, 
does save considerable time, and more so as t grows and the posterior for the 
parameter starts to change less. 
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uniform prior, without VIA 


FOFI 
POFI 


bias st. dev. MSE 
.1098 .6552 .4414 
.0424 .6376 .4083 


inaccurate prior, without VIA 


FOFI 
POFI 


bias st. dev. MSE 
.0123 .6458 .4172 
.0106 .6058 .3671 


fixed control (u t — .4 for all t) 


fixed 


.5336 .8298 0.9732 



Table 7: Simulation results for the 
POFI and FOFI and with and with 



uniform prior, with VIA 


FOFI 
POFI 


bias st. dev. MSE 
.0743 .5895 .3531 
.0145 .5970 .3567 


inaccurate prior, with VIA 


FOFI 
POFI 


bias st. dev. MSE 
.0054 .6035 .3643 
.0083 .6203 .3848 



/R Model using two kinds of priors, 
VIA. 



running time of VIA 
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" 9?&% u^i °o- 
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Time t 
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Figure 2: Running time of VIA at each time step t, for POFI using a uniform 
prior for the PCR model. 
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6 Discussion 



In this paper we compared two possible ways to conduct experimental design 
in parametric POMDP's, based on using dynamic programming to maximize 
either the Partial Observation Fisher Information or the Full Observation 
Fisher Information. Settings can arise where controls chosen by FOFI are not 
optimal, due to focusing on the underlying process rather than the observed 
process, and in these cases controls chosen with POFI often perform better, 
as in the six state example. In many of the examples analyzed they seemed 
to perform similarly. 

In recent years, there has been growing interest in statistical procedures 
within dynamical systems, such as parameter estimation and hypothesis test- 
ing, and many of these procedures could be performed more efficiently given 
good experimental design. In this paper we fully discretized the state and 
observational spaces to transform dynamical systems with stochastic errors 
into partially observed Markov decision processes, allowing us to use the 
methods developed for POMDP's to our advantage. 

We also noted how the problem of parameter dependence can be overcome 
by averaging over a prior. Additionally given that there is enough time 
between consecutive time steps, we showed how the controls can be efficiently 
updated online using observations gathered so far, by using a variant of the 
Value Iteration Algorithm. This was demonstrated in the PCR example. 

Finding controls that maximize information about parameters is a com- 
putationally challenging task. We have successfully demonstrated techniques 
for up to two dimensional systems, for a one dimensional parameter. Adding 
dimensions in state, parameter or observation space quickly make the meth- 
ods considered computationally intractable. Considering a longer lag of past 
observations for POFI might also increase accuracy, but again at the cost of 
computation time. In order to extend these methods to higher dimensions 
one could possibly use techniques of approximate dynamic programming, see 
Powell @]. 
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7 Appendix 



7.1 Calculations for the Dynamic program 

In order to run the POFI dynamic program or the VIA under POFI, we 
need to have p(y t+1 \y t , u t) y t -i, itt-i, 0) and logp(y t+ i|y t , u t , y t _ u u t ^, 0) 
calculated for every possible combination of y t ,u t ,y t -i,u t -i. Let K denote 
the size of the sample space of x t , L denote the size of our observation space, / 
the number of controls one can use. Thus we need to calculate 2L 3 / 2 numbers. 
We have that 

p(y t +i\yt,ut,y t -i,ut-i,9) 



= ^2^2p(y t +i\xt+i)p(xt+i\x t , u t , 0)p(x t \y t , yt-i, «t-i, 0) 

Xt+l x t 

and using Bayes rule: 

v(?t\Vt,Vt-\,ih-\,&) 
and again 



p{x t -i\y t -i 



P(yt\x t ) Ex t „! P(x t \x t -i, ut-i, 0)p(x t - 1 \y t - 1 ) 
Ex* P(Vt\ x t) Ex«_i p(x t \x t -i,ut-i, 0)p(x t -i\y t -i) 

p(y t -i\x t -i)p(x t -i) 



^PiVt-ilxt-Mxt-i) 



and we assume a uniform prior for x t _i. To calculate — logp(y t+ i\y t , u t , yt-i, Ut-i, 1 
we take the log of p(y t+ i\y t ,u t ,y t -i,u t -i,0) and take derivates by carrying 
out the chain rule through the sums above. The calculations are given more 
explicitly below in array form. 

We begin by constructing a L x K matrix E such that {Ej ti = p(y t = 
y j \x t = x 1 )} and a K x K diagonal matrix B j for each state in the observed 
y t state space. 



p{yt 



y 3 \x t 



x 1 ) 



p(yt 



y J \x t 



x K ) 



Let K 1 be a K x L matrix, such that K l = {K} • } = {p(x t = x H \y t = y^ 1 )}- 



We then let K 2 be a K x L x L x / array such that \K\ 



»2,J2Ji,ri 



} = {p(x t = 
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x l2 \yt = U^jUt-i = y^jUt-i — u Tl ,9) which can be calculated by fixing 
yt = yi 2 and u t -i = u ri and then 

yC^ — 

l:K,j 2 ,l:L,n \B^ 2 P ri K l 

Now let Q be a L x L x / x L x I array such that 

Qh,h,r2,h,n =p(yt+i = y J3 \yt = y n ,u t = u r2 ,y t -i = y n ,u t -i = u ri ,e) 

which can be calculated by fixing y t = y^ 2 ,u t = u r2 ,u t -i = u Tl and then 

Ql-L,j 2 ,r 2 ,j 3 ,ri — EP V2 Kl:K,j2,l:L, ri 

Now let C be a L x L x / x L x I array such that 
d 2 

c i3j a ,r a ji,n = -7^2 logp(j/ t +i = = y j2 ,^ = u T2 ,y t - X = y 3 \u t -i = u ri ,9) 

We calculate this by fixing y t+ \ = y-* 3 ,yt = y-* 2 ,u t = u r2 ,u t -i = u ri and 
setting = K^. K j 1:L which will be a K x L matrix, and then 

lB^{P r2 Kj + 2P r2 K 2 m + P r2 K 2 m )lB^P r2 K 2 m - {lB^{P r2 Kj + P r2 Kj)) 2 

We need the derivatives of K 2 or p(x t \yt, yt-i, Ut-i) and those can be calcu- 
lated by fixing y t = y^ 2 and u t -i = u Tl and then 

• 2 B^iP^K 1 + P^K^ilB^P^K 1 ) - B h P^K 1 {\B h {P T ^K x + F 1 ^ 1 )) 



(lB^PnK 1 ) 



2 



*™ = (lB^K^ {[Bh{priKl + + prikl )( 1BnpriKl ) 

_ B hpri K i( 1B h(pri K i + 2 p^ 1 + P ri K 1 ))](l J B i2 P ri K 1 ) 

^[P^P^+P^XlP^P^ 1 )^ 

To run the POFI dynamic program, or VIA under POFI, we can now use 
the arrays C and Q as described in the corresponding pseudocodes. 
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